Методи та алгоритми комп'ютерного рішення
диференціальних рівнянь
Введення
Для того, щоб описати динаміку різних процесів, що протікають у природних і в технічних системах, складають, спираючись на фізичні закони, диференціальні рівняння. Так, зокрема, доводиться робити при дослідженні функціонування автоматичних систем; роботи суднових енергетичних комплексів, електричних агрегатів, суднових допоміжних механізмів, систем навігації і т.д. У ряді випадків ці рівняння допускають лінеаризацію і можуть бути записані у вигляді:
де y (t) - невідома функція, a 0, a 1, ... a n - постійні коефіцієнти, а j (x) - деяка відома функція незалежного аргументу t, яка зазвичай висловлює зовнішній вплив, який чиниться на систему.
1. Мета контрольної роботи
Придбання навичок алгоритмізації та програмування задач чисельного інтегрування звичайних диференціальних рівнянь і систем з наступним моделюванням результатів на персональному комп'ютері і поданням їх у вигляді таблиць і графіків.
В результаті виконання контрольної роботи студент зобов'язаний:
1. Навчитися вирішувати лінійні диференціальні рівняння чисельними і символьними методами в рамках пакета комп'ютерної математики MathCAD.
2. Ознайомитися з основними алгоритмами існуючих комп'ютерних методів.
3. Визначити точність цих методів шляхом порівняння результатів, отриманих шляхом наближеного й аналітичного розв'язків.
2. Аналітичні методи
Загальне рішення диференціального рівняння n-го порядку - невідома функція y (t) - містить n довільних постійних. Їх можна визначити, знаючи початкові умови, що накладаються на невідому функцію і на її похідні до (n-1)-порядку включно. Аналітично (у символьному вигляді) такі рівняння вирішують класичним і операційним методами.
2 .1 Класичний метод
В обмеженому числі випадків виду лівій частині (1) допускає таке перетворення, яке дозволяє знайти рішення шляхом безпосереднього інтегрування, однак у загальному випадку порядок рішення - іншою.
Рішення неоднорідного диференціального рівняння (з ненульовою правою частиною) є сумою загального розв'язку відповідного однорідного диференціального рівняння y 1 (t) та приватного рішення y 2 (t) неоднорідного диференціального рівняння (1).
Рішення однорідного рівняння шукаємо у вигляді:
яке має n коренів -
де С i - довільні постійні, які знаходяться з початкових умов.
Є правила, що дозволяють визначити вид y 2 (t) приватного рішення в залежності від виду правої частини - функції j (t). Подальша підстановка спільного рішення у вихідне диференціальне рівняння дозволяє знайти невизначені константи C i у виразі для y 1 (t).
«Класичний» метод аналізу процесів в даний час використовується тільки у випадку найпростіших систем, оскільки необхідність знаходження приватного рішення часто призводить до складних перетворень, а також, крім рішення характеристичного рівняння додатково необхідно скласти і вирішити n рівнянь для визначення сталих інтегрування.
2 .2 Метод операційного числення
Суть методу полягає в проведенні інтегрального перетворення Лапласа функції, що входить до складу диференціального рівняння, за правилом:
де s = a + j × b - комплексна змінна величина.
Це перетворення зіставляє функції дійсного змінного функцію комплексного змінного
Перетворення Лапласа можна виконати, використовуючи блок символьних обчислень MathCAD. Цей же блок дозволяє виконати і зворотне перетворення Лапласа, відповідно до співвідношення:
де
Перетворення Лапласа зводить диференціювання функції оригіналу до множення її образу на комплексну змінну s, тому рішення диференціального рівняння в просторі оригіналів зводиться до вирішення алгебраїчного рівняння в просторі зображень.
Порядок вирішення диференціального рівняння з допомогою операційного числення представляється наступним:
- Виконуючи перетворення Лапласа лівої і правої частини диференціального рівняння, враховуємо початкові умови і переходимо від диференціального рівняння для функції оригіналу y (t) до алгебраическому рівнянню для Лаплас образу - Y (s);
- Вирішуючи рівняння алгебри, знаходимо в просторі зображень в явному вигляді вираз для Y (s);
- Виконуючи зворотне перетворення ми знаходимо невідому функцію y (t).
Всі етапи цієї процедури можуть бути автоматизовані і виконані в рамках пакету MathCAD (приклад 1).
Слід зауважити, що пакет MathCAD далеко не завжди здатний виконати в символьній формі результат зворотного Лаплас перетворення. Справа в тому, що в блок символьних перетворень пакету закладені правила виконання даної процедури для виразів записаних у вигляді елементарних дробів. Тому Лаплас образ попередньо розкладається на елементарні дроби. Однак, якщо корені полінома в знаменнику представляються у вигляді комбінації складних радикалів, то MathCAD «відмовляється» працювати. У цьому випадку йому необхідно «допомогти» врукопашну виконавши розкладання полінома в знаменнику відповідно до співвідношення:
де s 1, s 2, ... s n - корені рівняння
Розглянута методика знаходження аналітичного рішення диференційних рівнянь може бути поширена на завдання вирішення системи диференціальних рівнянь. У цьому випадку необхідно вирішити не одне алгебраїчне рівняння для Лаплас-образів, а систему алгебраїчних рівнянь за допомогою тієї ж процедури блоку рішень Given - Find. Відзначимо, що на відміну від систем комп'ютерної математики Mathematica 2.2.2 і Maple V R3/R4, які легко дозволяють аналітично вирішити лінійне диференціальне рівняння з допомогою вбудованих засобів. Система MathCAD передбачає «ручні процедури» запуску прямого перетворення Лапласа, складання за його результатами алгебраїчного рівняння і, після його рішення, запуску процедури зворотного перетворення Лапласа.
3. Чисельний метод рішення диференціальних рівнянь
Часто при аналізі вивчаються в технічних або в природних системах процесів доводиться враховувати наявність нелінійного поведінки функцій, що описують характеристики їх елементів. Це в свою чергу визначає появу нелінійностей в диференціальних рівняннях, які тепер вже не можуть бути записані у формі (1). Наявність нелінійностей в диференціальних рівняннях обумовлює неможливість їх точного аналітичного рішення, а наближені аналітичні методи часто призводять до громіздким викладкам.
Крім того, коефіцієнти в лівій частині диференціального рівняння можуть бути визначені з експерименту з помилкою, що, значною мірою, знецінює отримується точний аналітичний результат.
І нарешті, точні методи не придатні для випадку, якщо права частина диференціального рівняння представлена не в аналітичній формі, а у вигляді таблиці або графіка.
У всіх цих випадках прекрасно продовжують «працювати» методи чисельного рішення. На відміну від аналітичних, вони дозволяють отримувати шукані залежності для будь-якої з описаних вище ситуацій. Алгоритми існуючих методів чисельного рішення були розроблені порівняно давно, проте поштовх до їх застосування був обумовлений розвитком обчислювальної техніки. Кожен з існуючих чисельних методів припускає заміну похідної на кінцеве прирощення і перетворення диференціального рівняння в рівняння в кінцевих різницях.
З цією метою інтервал пошуку рішення розбивається на безліч відрізків і рішення шукається на кожному з цих шматочків. Ясно, що чим дрібніше крок розбиття, тим точніше виходить результат. Тому, ефективне застосування чисельних методів (при вирішенні реальних, а не навчальних задач!) Припускає використання комп'ютерів з достатнім швидкодією.
Використання для чисельного рішення диференціальних рівнянь комп'ютерного пакету MathCAD передбачає знання алгоритму роботи чисельних методів для розумного їх застосування (знання меж застосування, оцінки точності, витрат комп'ютерних ресурсів та ін.) Справа в тому, що до результатів комп'ютерних обчислень завжди потрібно ставитися критично; аналізувати їх на правдоподібність, і для того, щоб уникнути "підводних каменів" при використанні будь-якого стандартного пакета, що реалізує чисельні методи, потрібно мати хоча б мінімальне уявлення про те, який саме чисельний метод реалізований для вирішення тієї чи іншої задачі, як він «працює» .
Оскільки для оцінки точності рішення необхідний матеріал для порівняння пропонується розглядати наближені (чисельні) методи рішення стосовно до тих завдань, для яких раніше було отримано аналітичне рішення (тобто лінійних рівнянь виду (1)).
Позначимо
Тому, ідея застосування чисельних методів для вирішення рівняння старших порядків, в принципі, нічим не відрізняється від ідеї чисельного рішення рівняння першого порядку, які розглядалися в лабораторній роботі № 10.
Неважко розширити застосування описаної вище методики на випадок системи лінійних рівнянь. У прикладах 2 і 3 представлені реалізації методу Ейлера в рамках векторної процедури і з допомогою програми-функції. В обох випадках отримано однакові результати, які повіряються за отриманим раніше (див. приклад 1) аналітичного розв'язування. Рішення наближене і точне сильно відрізняються і величина відносної помилки (для обраного кроку) складає ~ 13%. Приклад 3 неважко оптимізувати і надати йому більш компактний вигляд, вважаючи початкові умови і праві частини системи рівнянь компонентами деяких векторів. Також припустима доопрацювання програми-функції на випадок будь-якого числа рівнянь шляхом зміни числа аргументів програми-функції, а також числа рядків у програмі: в завданні початкових умов і в циклі обчислення масиву рішень.
Аналогічно методом Ейлера векторна і програмна реалізації обчислень за методом Рунге-Кутта можуть бути поширені на випадок рішення системи диференціальних рівнянь. Пряме рішення завдання в цьому випадку є досить громіздким і бажано передбачити розрахунку констант методу в рамках внутрішніх процедур. У пакеті MathCAD є вбудовані функції, вирішальні подібні завдання. Так, процедуру розрахунку наближеного рішення за методом Рунге-Кутта вирішується за допомогою функцій rkfixed і Rkadapt. Приклад вирішення тієї ж системи диференціальних рівнянь дано в прикладі 4. Відмінність в точності отриманих рішень для даного диференціального рівняння невелика, однак, якщо рішенням диференціального рівняння є сильно осцилююча функція, то метод зі змінним кроком забезпечує більшу точність.
Приклади рішення різних рівнянь за допомогою комп'ютерного пакету MathCAD.
| Приклад № 1 |
| Приклад № 2 |
| Приклад № 3 |
| Приклад № 4 |
Замінюючи функцію y (x) на Y, одержимо вираз, який потім спрощуємо за допомогою функції collect.
Отримуємо вихідне рівняння алгебри яке вирішуємо в символьному вигляді:
Проводимо перевірку правильності початкових умов, використовуючи функцію subtitute